'r  —  ” 

AD-A115  483  AIR  FORCE  INST  OF  TECH  WR16HT-PATTERS0N  AFB  OH  SCHOO— ETC  F/6  12/1 
smooth  fit  analysis  OF  MOESSBAUER  SPECTRA. (U) 
mar  82  T  U  BAZZOLl 

UNCLASSIFIED  AFIT/6NE/PH/82-2  NL 


I 


DEPARTMENT  OF  THE  AIR  FORCE 
ANt  UNIVERSITY  (ATC) 


ELECTE 

JUN14ia8Et 


AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


AFIT/GNE/PH/82-2 


SMOOTH  FIT  ANALYSIS  OF 
MOSSBAUER  SPECTRA 

THESIS 


AFIT/GNE/PH/82-2  Thomas  L.  Bazzoli 

Lieutenant  USAF 


ll 


\ 

V 


Approved  for  public  release;  distribution  unlimited. 


AFIT/GNE/PH/82-2 


SMOOTH  FIT  ANALYSIS  OF 
MOSSBAUER  SPECTRA 


THESIS 

Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 

in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science 


by 

Thomas  L.  Bazzoli,  B.S.  Chemistry 
Lieutenant  USAF 

Graduate  Nuclear  Effects 
March  1982 

Approved  for  public  release;  distribution  unlimited. 


Acknowledgments 


I  would  like  to  thank  and  acknowledge  my  advisor, 

Dr.  Donn  Shankland  for  suggesting  this  project  and  for 
pointing  me  in  the  right  direction  when  I  had  wandered  off 
on  my  own  tangents.  I  would  also  like  to  extend  special 
thanks  to  my  classmates  of  GNE-82M  who  were  eager  to  encour¬ 
age  me  when  my  problems  seemed  insurmountable. 


ii 


Contents 

Page 

Acknowledgments  .  ii 

List  of  Figures .  v 

List  of  Tables . vii 

Abstract . viil 

I.  Introduction  .  1 

Background  .  1 

Method  .  4 

Purpose  .  7 

Overview  .  7 

II.  Development  of  Fitting  Technique  .  8 

Least-Squares  Approach  .  8 

Limitations  on  the  Magnetic  Field 

Density . 10 

Nonnegative . 11 

Smooth . 11 

Minimization  .  12 

Numerical  Solution  for  the  Integral  .  16 

Lineshape  Function  .  16 

Set  of  Functions  9>n . 18 

Integral  Approximation  .  20 

Recursion  Relation  .  24 

III.  Computer  Program  Development  .  30 

Computation  of  rQ(c)  .  30 

Subroutine  FIX . 32 

Subroutine  CONS . 33 

Subroutine  INTEG  .  35 

Subroutine  RITEO  .  36 

Determination  of  p(H) . 36 

Subroutine  READN  .  36 

Subroutine  MATRl  .  39 

Subroutine  HERMIT  .  39 

Subroutine  MESCHP  .  39 

Subroutine  MAGFD  .  40 

Subroutine  PLTMAG  .  40 

Computation  of  Fit  Spectrum . 41 

Subroutine  SPCTRM  .  41 

Subroutine  PLTSPT  .  42 


iil 


Page 


Other  Routines . 42 

Subroutine  CONST  .  42 

Subroutine  MATRN  .  43 

IV.  Program  Validation . 44 

V.  Results,  Recommendations  and 

Conclusions . 59 

Optimizing  Results  .  59 

Conclusions . 60 

Recommendations  .  62 

Bibliography . • . 64 

Appendix  A;  Input  Parameters  for  Program 

RPARAM . 66 

Appendix  B:  Computer  Program  Listings  .  70 

Vita . 99 


iv 


List  of  Figures 

Figure  Page 

1  Magnetic  Field  Density  and  the  Resulting 

Spectra  from  an  Amorphous  Material  .  2 

2  Magnetic  Field  Density  and  the  Resulting 

Spectra  from  an  Alloy  in  the 

Crystalline  State  .  3 

3  Flow  Diagram  for  RPARAM . 31 

4  Generalized  Flowchart  for  MOSFIT  pt.  I  .  .  ,  .  37 

5  MOSFIT  Flowchart  pt.  II . 38 

6  Examination  of  the  General  Lineshape 

Produced  With  Given  Data . 45 

7  Magnetic  Field  Density  Used  to  Generate 

Test  Spectrum . 49 

8  Input  Test  Spectrum . 50 

9  Magnetic  Field  Density  from  Initial 

Test  Run . 52 

10  Fit  Spectrum  from  Numerical  Integration 

in  Initial  Run . 53 

11  Magnetic  Field  Density  With  2x 

Correction  Factor  .  54 

12  Fit  Spectrum  by  Numerical  Integration 

and  2x  Factor . 55 

13  Fit  Spectrum  Using  r^Cc)  Values  Without 

2x  Factor . 56 

14  Fit  Spectrum  Using  rjj(c)  Values  and 

2x  Factor . 57 

15  Result  Optimization  .  60 

A-1  Relation  of  Input  Parameters  to 

Spectrum  . . 66 

A-2  Splitting  of  the  Fe^"^  Nucleus . 68 


V 


Figure 


B-1 

Flowchart 

for 

Subroutine 

CONS  . 

B-2 

Flowchart 

for 

Subroutine 

INTEG  . 

B-3 

Flowchart 

for 

Subroutine 

COEFF  . 

B-4 

Continuation  of  COEFF  Flowchart 

B-5 

Flowchart 

for 

Subroutine 

HERMIT 

B-6 

Flowchart 

for 

Subroutine 

MAGFD  . 

B-7 

Flowchart 

for 

Subroutine 

SPCTRM 

vl 


List  of  Tables 


Table  Page 

57 

I  Input  Parameters  for  Fe . 32 

II  Evaluation  of  FNJ . 46 


vii 


AFIT/GNE/PH/82-2 


Abstract 

A  method  of  fitting  the  magnetic  field  distribution 
from  Mossbauer  spectra  using  a  linear  least-squares  tech¬ 
nique  was  proposed  where  the  distribution  was  assvimed  to  be 
an  expansion  in  Hermite  polynomials.  This  method  was 
developed,  programmed  and  tested  on  a  test  spectriom  with  a 
known  distribution  of  magnetic  field.  The  programming  was 
not  entirely  successful  in  reproducing  this  test  spectra 
and  no  actual  data  was  analyzed  to  support  the  effectiveness 
of  this  method. 


vili 


SMOOTH  FIT  ANALYSIS  OF 
MOSSBAUER  SPECTRA 

I .  Introduction 

Mossbauer  spectroscopy  is  known  for  its  applicability 
in  many  fields  and,  as  a  result,  has  become  a  valuable  tool 
in  most  disciplines  of  science.  Today,  the  principles 
involved  in  nuclear  gamma-ray  resonance  spectroscopy  are 
well  documented  in  the  literature  (Ref  1,2).  However,  while 
a  Mossbauer  spectrum  can  provide  a  wealth  of  information  on 
the  properties  of  a  material  under  investigation,  much  of 
this  information  can  be  lost  in  the  statistical  fluctuations 
of  the  data  and  the  complexity  of  the  spectrum.  In  response 
to  the  problems  of  Mossbauer  analysis,  many  researchers  have 
developed  techniques  of  fitting  Mossbauer  spectra.  This 
work  examines  a  fitting  procedure  developed  by  Dr.  Donn  G. 
Shankland  (Ref  3). 

Background 

A  Mossbauer  spectrum  can  provide  a  means  to  investi¬ 
gate  the  microscopic  fields  of  a  material .  One  class  of 
materials  that  is  of  interest  to  the  Air  Force  is  amorphous 
metal  alloys,  which  are  known  for  their  flexibility  and 
high  strength.  The  disorganized  lattice  structure  of  these 
alloys  results  in  various  magnetic  fields  from  short  range 
ordering  (Ref  4).  These  fields  are  distributed  in  a 
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probability  distribution  commonly  called  the  magnetic  field 
density  and  represented  as  pCH)  .  In  the  amorthous  state, 
this  distribution  is  continuous  and  the  resulting  Mossbauer 
spectrum  is  broad  without  distinct  peaks.  An  example  of 
the  distribution  and  spectrum  expected  from  an  amorphous 
material  is  shown  in  Figure  1. 


Figure  1.  Magnetic  Field  Density  and  the  Resulting 
Spectra  from  an  Amorphous  Material 

On  the  other  hand,  as  the  amorphous  alloy  crystallizes, 

distinct  magnetic  domains  will  develop — each  with  a  discrete 

magnetic  field.  The  alloys  of  interest  to  the  Air  Force 

57 

consist  of  various  iron  compounds.  When  Fe  is  present  in 
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an  absorber  and  a  magnetic  field  exists,  the  resulting 
Mossbauer  spectrum  will  consist  of  six  peaks  that  correspond 
to  the  six  allowed  transitions  in  this  isotope.  Each  dis¬ 
crete  magnetic  field  will  cause  this  spread  of  the  Mossbauer 
lines.  As  a  result,  the  Mossbauer  spectrum  of  the  crystal¬ 
line  state  will  consist  of  a  set  of  overlapping  six  peak 
spectra,  each  set  corresponding  to  a  discrete  magnetic  field. 
A  conceptual  representation  of  the  magnetic  field  density 
and  spectrum  of  the  crystalline  alloy  is  found  in  Figure  2. 


Figure  2.  Magnetic  Field  Density  and  the  Resulting 
Spectra  from  an  Alloy  in  the  Crystalline 
State 
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This  study  investigates  a  method  for  finding  the 
magnetic  field  density  from  amorphous  metal  alloys  and  using 
the  result  in  a  fit  of  the  Mossbauer  spectrum. 

Method 

The  data  collected  in  the  investigation  of  the 
Mossbauer  effect  is  generally  in  the  form  of  a  digitized 
spectrum  that  represents  the  interaction  of  gamma  radiation 
with  a  given  detector.  Since  random  processes  are  involved 
in  radioactive  decay,  absorption,  reemission  and  detection, 
the  individual  data  points  are  subject  to  a  degree  of  sta¬ 
tistical  fluctuation.  The  predominant  method  of  minimizing 
this  statistical  error  involves  the  application  of  a  least- 
squares  fitting  technique. 

In  general,  least-squares  fitting  requires  the  compar' 
ison  of  the  data  with  a  functional  approximation  of  this 
data.  The  approximating  function  should  be  chosen  so  that 
the  sum  of  the  squares  of  the  difference  between  the  data 
and  this  function  is  minimized.  That  is,  if  the  data  is 
represented  as  n^  and  the  functional  approximation  is  nCc), 
then  the  expression 

S  =  ZCn^  -  n(c))^ 

is  minimized.  This  minimization  is  accomplished  by  deter¬ 
mining  the  parameters  in  nCc)  which  can  be  varied,  and  then, 
the  partial  derivatives  of  s  with  respect  to  each  of  these 
parameters  are  found  and  set  equal  to  zero. 
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The  preliminary  task  in  performing  the  least-squares 
fitting  is  to  develop  the  approximating  function  or  theoret¬ 
ical  model  to  describe  the  data  obtained.  This  model  is 
used  to  represent  the  expected  or  "true"  values  to  be  com¬ 
pared  with  the  experimental  data  in  the  analysis.  When  a 
multichannel  analyzer  is  used  to  collect  a  Mossbauer  spec¬ 
trum,  then  the  number  of  counts  in  channel  c,  nCc),  is 
found  by 


n(c)  =  ngCc)  -  ^/‘“dHp(H)L(c.H) 

where  HgCc)  is  the  number  of  background  counts  in  channel  c, 
H  is  the  magnetic  field  strength,  p(H)  is  the  aforementioned 
magnetic  field  density,  and  L(c,H)  is  a  function  represent¬ 
ing  the  lineshape  of  the  spectral  peaks  (Ref  4). 

However,  problems  tend  to  arise  in  the  evaluation  of 
this  function.  The  approximating  equation  is  a  Fredholm 
integral  function  of  the  first  kind  and  problems  occur  when 
the  integral  is  approximated  numerically  (ilef  5).  The 
solution  of  this  integral  is  nonunique  and  the  problem 
becomes  ill-posed  since  the  solution  process  is  highly 
sensitive  to  the  data.  As  a  result,  it  is  advantageous  to 
find  an  analytical  solution  to  this  problem  by  choosing  the 
appropriate  functional  forms  of  L(c,H)  and  pCH). 

K7 

Since  Fe  is  the  absorber  in  the  amorphous  metal 
alloys,  it  can  be  assumed  that  the  spectrum  of  the  alloy 
will  be  fundamentally  related  to  the  six  peak  spectrum  of 
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Fe’’  in  the  presence  of  a  magnetic  field.  Therefore,  the 
overall  lineshape  function  can  be  represented  as  the  sum  of 
six  individual  peak  shapes.  Most  work  in  Mossbauer  spec¬ 
troscopy  describes  individual  absorption  peaks  with  a  simpl 
Lorentzian  formula.  In  contrast,  this  work  describes  the 
lineshape  as  a  sum  of  six  Gaussian  peaks  in  order  to  facili 
tate  the  analytical  solution  of  the  integral. 

In  the  simplification  of  the  integral,  other  inves¬ 
tigators  have  assumed  a  shape,  or  modelled  the  magnetic 
field  density  in  their  attempts  to  fit  Mossbauer  spectra 
(Ref  6-10).  Shankland  chose  to  fit  p(H)  by  an  expansion  in 
Hermite  polynomials.  Combining  this  fit  of  the  magnetic 
field  density  with  the  lineshape  function  reduces  the  inte¬ 
gral  can  be  approximated  by  a  given  number  of  series  terms. 
The  solution  is  developed  so  that  the  number  of  series 
terms  chosen  will  correspond  to  the  number  of  Hermites  in 
the  expansion  of  p(H).  The  integral  terms  are  incorporated 
in  a  matrix  minimization  problem  in  the  least-squares 
analysis.  The  solution  of  the  minimization  problem  will 
yield  a  magnetic  field  distribution  that  is  continuous  and 
best  represents  the  microscopic  field  distribution  that 
produced  the  data. 

Two  problems  may  arise  in  this  approach.  If  a  sharp 
peak  (or  peaks)  occurs  in  the  magnetic  field  density,  a 
large  number  of  polynomials  will  be  required  to  reproduce 
the  peak.  Therefore,  errors  can  be  introduced  when  only  a 
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few  Hermites  are  used  in  the  expansion  of  p(H).  In  addi¬ 
tion,  the  p(H)  generated  may  consist  of  nonphysical  negative 
values.  This  requires  the  addition  of  the  constraint  that 
the  magnetic  field  always  be  positive.  This  study  attempts 
to  fit  Mossbauer  spectra  by  representing  p(H)  as  an  expan¬ 
sion  of  Hermite  polynomials. 

Purpose 

The  purpose  of  this  study  was  to  investigate  the 
effectiveness  of  using  Dr.  Shankland's  method  for  fitting 
p(H)  and  subsequently  using  the  result  to  determine  the  fit 
spectrum.  Implementation  of  this  method  involved  writing 
two  computer  codes  using  the  numerical  procedure  developed, 
validating  the  method  by  generating  a  test  spectrum  with 
known  p(H)  and  then  applying  the  method  to  actual  Mossbauer 
data. 

Overview 

The  theory  of  Mossbauer  spectroscopy  is  well  known 
and  will  not  be  discussed  (Ref  1,2).  Chapter  II  contains 
a  detailed  development  of  Dr.  Shankland’s  fitting  procedure. 
The  implementation  of  this  theory  is  discussed  with  the 
development  of  the  computer  codes  in  Chapter  III.  Chapter 
IV  describes  the  procedures  used  in  the  validation  and  test¬ 
ing  of  the  programs.  The  final  chapter  includes  the  results, 
conclusions  and  recommendations. 


7 


II .  Development  of  Fitting  Technique 

In  the  measurement  of  the  Mossbauer  effect,  a  digital 
spectrum  consisting  of  background  counts  and  a  number  of 
absorption  peaks  is  produced.  This  digitized  data  results 
from  gamma  ray  emission,  absorption  and  detection.  Conse¬ 
quently,  these  measurements,  or  counts  in  the  experimental 
sense,  are  subject  to  inherent  statistical  fluctuation  that 
introduces  uncertainty  into  the  results.  In  order  to 
investigate  the  microscopic  fields  of  amorphous  metal  alloys, 
the  magnetic  field  distribution  that  produces  this  data  must 
be  found.  Dr.  Shankland  has  developed  a  method  of  fitting 
the  magnetic  field  distribution  of  Mossbauer  spectra  by 
employing  a  linear  least-squares  technique. 

Least-Squares  Approach 

Since  the  experimental  data,  n^i^,  is  subject  to  the 
errors  of  counting  statistics,  it  is  desirable  to  find  a 
"best  fit"  of  the  data  with  respect  to  the  "true"  or 
expected  results.  A  common  means  of  curve  fitting  is  the 
method  of  least-squares  (Ref  11).  In  general,  a  least- 
squares  fit  chooses  certain  parameters  to  minimize  the 
relation 

I  w.v  2  (2.1) 

1  ^  ^ 
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where 


represents  a  weighting  factor 

is  the  residual  or  difference  between  the  true 
and  observed  results 

When  a  multichannel  analyzer  is  used  to  collect  the  Mossbauer 
spectra,  the  "true"  results,  n(c),  can  be  predicted  by  the 
approximating  function 

n(c)  =  ngCc)  -  Q/“dHp(H)L(c,H)  (2.2) 

where 

ng(c)  is  the  number  of  background  counts  in 
channel  c 

pCH)  is  the  magnetic  field  density  from  all  sources 
LCc,H)  is  a  lineshape  function  used  to  describe  the 
absorption  peaks  in  the  spectrum 
With  the  residual,  =  n(c)  -  ,  Eq  (2.1)  becomes 

Z  w. (n(c)  -  n. )^ 

i  ^  ^  (2.3a) 


or 

Z  w.[n.  -  ngCc.)  +  ^/*dHp(H)L(c^,H)]2  (2.3b) 

In  counting  statistics,  errors  are  found  to  vary  like 
©i  =  /n^^  (Ref  12).  If  the  error  is  examined  as  a  fraction 
of  the  Individual  data,  /n^/n^^  ,  it  is  evident  that  this 

fractional  error  can  be  found  to  vary  as  ^ 


9 
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Since  is  a  measure  of  the  expected  deviation  from 
the  true  value  for  a  given  data  point,  it  would  be  reason¬ 
able  to  minimize  the  residual  weighted  by  — ^  (Ref  11). 
Hence,  the  least  squares  relation  is 

So  =  I  r  n^[n^  -  ng(c)  +  Q/"dHp(H)L(c ,H) ] ^  (2.4) 

where  the  ^  is  a  constant  of  proportionality.  Finally,  the 
background  ng(c)  must  be  described.  Since  the  background 
is  diffuse  and  is  assumed  to  be  curved,  it  can  be  adequately 
modelled  by  a  quadratic  of  the  form 

ng(c)  *  +  b^c  +  bgC^  (2.5a) 

or 

2  Y 

nn(c)  =  Z  be'  (2.5b) 

®  Y=o  ^ 

Thus,  Eq  (2.4)  becomes 

=  -J  ^  ^  ,^/“dHp(H)L(c,H)]2  (2.6) 

o  y  ^ 

which  is  the  system  to  be  minimized. 

Limitations  on  the  Magnetic  Field  Density 

It  is  also  desired  that  this  fit  be  subject  to  the 
criteria  that  the  magentic  field  density,  p(H),  be  smooth 
and  nonnegative.  The  fit  must  be  subject  to  these  condi¬ 
tions  because  the  smoothness  of  the  fit  spectra  is  directly 
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related  to  the  smoothness  of  p(H).  In  addition,  it  is 
important  to  eliminate  nonphysical  negative  values  of  the 
distribution. 


Nonnegative .  Generally,  there  may  be  one  or  more  of 
the  discrete  points  where  p(H^)  will  solve  as  negative. 
These  points  must  be  constrained  to  zero  to  fulfill  the 
condition  that  p(H),  which  represents  a  probability  distri¬ 
bution  of  magnetic  field  strengths,  be  positive  or  zero. 
With  the  addition  of  constraints 


S_  = 


I  I  “i'“i 


S  b  +  /"dHpCH)L(c,H)l^ 
y  y 


-  2  >  (2.7) 

q  q 


where 

p(H^  )  are  the  negative  values  of  p(H^) 

q 

Xq  are  the  Lagrangian  multipliers  for  each 

constrained  point 

Smooth .  The  roughness  of  pCH)  could  be  measured  by 
a  quadratic  form  in  p(H)  and  its  various  derivatives. 
Therefore,  a  smoothing  function  can  be  represented  as 


^o  =  o-'’"<i”P(H)g(^)pCH) 


C2.8) 


-d 

with  g('^— *■)  as  an  operator  in  the  derivatives  of  p(H). 
dH^ 
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However,  the  appropriate  operator  is  dependent  on  the  func¬ 
tional  form  of  the  magnetic  field  density.  As  a  result,  it 
is  desirable  to  represent  the  smoothness  function  as 

Ro  =  o/“dHp(H)g(-A2)p(H)  (2.9) 

2 

where  the  g(-A  )  represents  any  operator  that  will  provide 

a  measure  of  roughness  for  the  specific  form  of  p(H). 

The  Rq  term  in  Eq  (2.9)  is  weighted  and  added  to  the 

system  so  that  it  may  be  minimized  simultaneously  with  the 

constrained  data  error  S  .  That  is, 

c 

S  =  <rs  +  BR  (2.10) 

c  o 

If  the  weight  terms  are  chosen  such  that  «  =  1  and 

A 

6  «  ^  ,  the  system  to  be  minimized  becomes 

S  =  -I  E  n^[ni  -  r  b^c^  +  ^/*dHp(H)L(c.  ,H)]^ 

-  E  A  P(H„  )  +  I  ^/"’dHp(H)g(-A2)p(H)  C2.ll) 
q  ^  q 


Minimization 

In  order  to  compute  the  "best"  smooth-fit  of  the 
experimental  data,  Eq  (2.11)  must  be  minimized  with  respect 
to  the  parameters  p(H),  b^,  and  A^.  Thus,  on  differentia¬ 
tion. 
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3S 


=  0  =  E  n^In.  -  E  b^c’  +  dH^p(H' )L(c .  ,  H' )  L(c.  .H) 


-  Xqd(H-Hq)  +  Ag(-A^)p(H)  (2.12a) 


3S 


3p(H) 


0  =  E  n.[n.  -  E  +  ^/*dH'p(H' )L(c . . H' ) ] (-c^) 


(2.12b) 


3S 


31. 


0  =  -p(Hq) 


(2.12c) 


where  5(H-Hq)  in  Eq  (2.12a)  represents  a  Dirac  delta  func¬ 
tion  at  Hq .  With  rearrangement,  the  three  equations  in 
Eq  (2.12)  yield 


/dH'[E  n^L(c.  ,H)L(c.  ,H')lp(H-')  +  Ag(-A"')  (H") 


-  E  [E  n.L(c^,H)cY]b^  -  lq5(H-Hq) 


-  E  n^^  L(c^,H)  =  0  (2.13a) 


fdH'lZ  n^c^L(c.  ,H')]p(H')  +  E  [E  n-cj'^^  ]b 
i  -y'  i  ^ 

-  E  n.^cT  =  0 
i  ^  ^ 


(2.13b) 


-  p(Hq)  =  0 


(2.13c) 


If  it  is  assumed  that  the  magnetic  field  density 
p(H)  can  be  expanded  in  a  complete  orthonormal  set  of 
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functions  0j,(H)  such  that 

p(H)  =  Z  a^p^(H)  (2.14) 

n 


and 


then  the  set  of  equations  in  Eq  (2.13)  becomes 

Z  n^L(c^,H)L(Cj.H')]  +  Ag(-A2)}^j^.(H-)a^- 

n  i 

-  1  6(H-H  )  -  E  [Z  n.L(c.  ,H)cpb 
^  ^  Y  i  ' 

-  Z  n^\(c^,H)  =  0  (2.16a) 

-Z^^/“’dH'[Z  n.c^L(c^,H")]pjj^(H')a^. 

Ti'  i 

+  Z  [Z  n.c -  Z  =  0  (2.16b) 

Yiii  Y 

'  ^,^n^”q^^n'  '  °  (2.16c) 


Notice  that  when  Eq  (2.16a)  is  multiplied  by  Pjj(H)  and  then 
integrated  over  H,  the  resulting  set  of  equations  can  be 
written  as 


Z  . 

_ ..  nn  n 
n 


-  Z  -  Z  C\„ 

y'  ny  Y  q  aQ  q 


-P. 


(2.17a) 
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(2.17b) 


-ZD..a^  +  EE  >.b..  =  q 
n  Y  n  YY  Y 

=0  (3.170 

Or,  in  matrix  form,  the  least-squares  minimization  problem 
is 

B  -D 

-D  E 

-C  0 

with 

®nn'  =  ;/dHdH'(8jj(H)[rnj^L(c^,H)L(c^,H') 


+  Ag(-A2)]P^.(H")  (2.19a) 

Cnq  =  (2.19b) 

=  /dH9)jjCH)[In.L(c.  ,H)c|]  (2.19c) 

YiY'  =  0,1,2  (2.19d) 

Pn  =  /dH(Z)jj(H)[In.2L(c^,H)3  (2.19e) 

=  En  2c  Y  y  =  0,1,2  (2.19f) 

Y  i  1  1 


— 

- 

a 

-P 

b 

= 

q 

X 

0 

MM  MM 

(2.18) 
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The  focus  of  the  problem  now  lies  in  finding  a  solu¬ 
tion  for  the  integral  I  =  /dH^5^(H)L(c ,H)  found  in  Eqs 
(2.19a),  (2.19c)  and  (2.19e). 

Numerical  Solution  for  the  Integral 

In  order  for  the  integral,  I,  to  be  solved  analyti¬ 
cally,  the  appropriate  forms  of  the  functions  p(H)  and 
L(c,H)  must  be  known. 

Lineshape  Function .  The  lineshape  function  L(c,H) 
consists  of  a  contribution  from  each  peak  that  occurs  in  a 
given  spectrum  such  that 

m 

L(c.,H)  =  r  f.if^(c.  ,H)  (2.20) 

j^l  ^  ^ 

where 

ra  =  number  of  peaks  in  the  spectrum 
f  .  =  recoil-free  fraction  or  intensity  of  the 

*J 

peak 

£.(c,H)  =  function  describing  the  actual  shape  of 
J 

individual  peaks 

In  general,  the  function  il.(c,H)  is  usually  assumed 
to  be  either  a  Lorentzian  or  Gaussian  lineshape  located  at 

Cj(H)  =  d^  +  d^  H  C2.21) 

where 

d^  is  the  channel  number  of  the  spectrum  center 
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is  a  factor  that  describes  the  "shift”  of  an 
individual  peak  with  magnetic  field 
c  .  is  the  resulting  channel  location  of  a  peak 

Therefore,  it  follows  that  when  H=0  ,  each  peak  will  be 

centered  at  d^ .  However,  the  peaks  will  spread  linearly, 
but  differently  with  varying  H. 

When  a  magnetic  field  causes  splitting  of  the  nuclear 
states,  the  Mossbauer  nucleus,  Fe®’,  has  six  allowed  transi¬ 
tions  which  will  yield  a  spectrum  consisting  of  six  peaks. 
Assuming  a  Gaussian  lineshape 


e  -  (2-h)^ 

^  2 

£.(H)  =  - -  j  =  1  to  6  (2.22) 

^  r  /2Tr 

with 

h  =  location  of  the  Gaussian  peak  center 
r  =  peak  linewidth 

By  substituting  Eq  (2.21)  for  h,  Eq  (2.22)  becomes 


Aj(c,do,d^,r,H) 


-(c-dQ-d^H)2 

_ 2r^ 

r 


(2.23) 


and  the  Integral  I  is 


I 


e  - - - 

Z  f  ./dH  - — -  P„(H) 

j  r  ° 


(2.24) 


Set  of  Functions  0  .  The  set  of  functions  0  must 
-  —  -  i-n  n 

combine  to  form  a  distribution  of  magnetic  fields.  Based 
on  Eq  (2.21),  it  is  clear  that  a  given  absorption  peak  can 
be  moved  with  varying  H.  In  general,  only  positive  field 
strength  is  considered,  but  by  allowing  a  negative  magnetic 
field  H  and  thereby  assuming  that  p(H)  is  symmetric,  an 
appropriate  choice  for  would  be  an  orthonormal  form  of 
the  Hermite  polynomials. 

Starting  with  the  generating  function  for  the  Hermite 
polynomials, 


^+2sZ  ^  ^  H^(2) 


(2.25) 


the  addition  of  a  weighting  factor  for  orthogonality, 
e”^  will  yield 


-s 


^+2sZ 


(2.26) 


The  set  of  functions  are  normalized  with  the  constant 


. —  — - —  such  that 

/2**n !  Vir 
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2  -s^/2+sZ/2" 

l/F" 


o  v/nT  /2Hn ! 


=  Z  ^  h(Z)  (2.27) 
o  /iTT  “ 


where  the  orthonormal  functions  0„(Z)  are 

n 

=  1 - = 

/tt~  /2"n!  (2.28a) 

with  the  generating  function 

^-Z^/2^-s^/2+sZ/5 

VF*  (2.28b) 


If  the  magnetic  field  density  is  redefined  in  terms  of  a 
dimensionless  variable,  Z,  then  p(Z)  =  Ia^0jj(Z)  with 
Z  =  H/Hq  .  Then  can  be  replaced  with  the  left  hand 

side  of  Eq  (2.27),  and  the  integral  becomes 


I 


-(c-dQ-d^H)2 

«  e  ^  -s2/2+sZ/2-z2/2 

E  f  .  /  dH  ^ ^ ^ - j - 

J  ^  (2.29a) 


which  is  related  to  the  functions  h^^  and  as 


»  ^n 


I  =  Z  f  Z  ^  dH  h  (Z)£.(c.d  di.r,H) 


(2.29b) 
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Integral  Approximation.  In  order  to  find  a  solution 
to  the  integral  in  Eq  (,2.29a),  the  exponential  terms  must 
be  combined  and  reduced.  First,  using  the  dimensionless 
variable  Z  defined  earlier,  where  H=H^Z  ,  the  integral  is 

I  =  E  f  j/H^dZ  e _ 2r^ _  .  J2+sZ/2-Z^  12 

^  T  /2Tr  I/tT  (2.30) 


Combining  the  exponential  terms  yields 


exponent  = 


„  fd.^H  ^  1  f  d.H  (c-d  )'| 


-  i  s2  + 

9 


(2.31) 


which  further  simplifies  as 

Z^-2Z[r^s/2-d^H^(c-d^^)] 

+  [T^s^+ic-d^)^] 

+  (d^H^)2 


Finally,  by  completing  the  square  the  exponent  is 


d?Ho^r= 


z  - 


r^s^^-djH^Cc-d^) 


L  J 


r^s^+Cc-d  )2 


r^^-diHoCc-dp)] 
r2  ^  (d^Hp)2 


(2.33) 


Since  the  last  two  terms  of  the  exponent  are  independent 
of  Z,  integration  of  the  exponential  yields  an  I  such  that 


Hp  exp 


I  =  Z  f  . 


j  ^  /r^+(d^Hp)^ 


1 

2 


+  (c-d 

o 


(r^s/2-d^Hp(c-dp))2 

~r^(r^+(d^Hp)^) 


I 


(2.34) 


Where  both  dj^  and  T  may  be  subscripted  j  to  correspond  with 
each  peak.  Assuming  that  the  linewidth  does  not  vary, 

Eq  (2.34)  can  be  rewritten  as 


H. 


I  =  r  f. 


^  /r^+(djHp)' 


exp 


(a^s^+bS+c')V 


(2.35) 


with 


1  - 


2r 


r^+(djHo)‘ 


(d^Hp)^-r2 

(d^Hp)^+r2 


(2.36a) 
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2/2  d^H^Cc-d^) 

^  ~  2  ^  2 

r  +(d^H^)  (2.36b) 

e-  .  .  (djH„)^(c-d^)^  _  (d-d„)^ 

r2  r^+(djH^)^  (2.36c) 


Therefore,  the  exponential  term  in  c'  can  be  placed  with 
the  constant  to  give 


H„e 


I  = 


°  2(r^+(dJH  )^)  ^  2  2 

E  f  .  -  - —  exp(--|  [a^s^+bs]  ) 

^  ^/iF  /f^+(d^H^^  (2.37) 


This  expression  is  then  defined  as 


(2.38a) 


which  is  equivalent  to  the  relation 


2  f.  E  —  „/”dH  h„(^)f  .(c,d^,d^,r,H)  (2.38b) 

-i^n/nT'  nH„j  ol 


so  that  the  number  of  functions  J  used  to  approximate  the 
integral  is  directly  related  to  the  number  of  Hermite  poly¬ 
nomials  in  the  description  of  p(H). 

If  the  exponent  in  Eq  (2.37)  is  expanded  in  a  series 

where 
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Recursion  Relation 


In  order  for  the  functions  to  be  generated  in  a 
computer  solution,  it  would  be  advantageous  to  find  these 
functions  by  recursion.  Since  the  Hermite  functions  satisfy 
the  recursion  relation 

2  X  Hj^CX)  =  H^+i(X)  +  2n  (2.42) 

and  recalling  that  the  Hermites  are  generated  by  the  rela¬ 
tion 


-s 


'+2s"X 


eo 

=  z 
o 


n ! 


Hn(X) 


(2.26) 


it  is  evident  that  the  functions  can  be  found  by  recur¬ 
sion 


-4(a^s^+bs)  _  -s'^+2s'X 

6  a  —  © 


(2.43) 


Solving  for  the  term  X  in  Eq  (2.43)  yields 


/(d^H^)^-r‘^  (2.44) 

with  a  and  b  as  defined  in  Eq  (2.36).  As  a  result,  from 
Eq  (2.26)  it  follows  that 


-i(a2s2+bs) 

e 


^-s'^-h2s'X  ^  j. 

o  /  2 


(2.45) 
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Furthermore,  if 


1 
2 


N(c.dJ,H^,d^,r)  = 


_  H^e  ^  r^+Cd^H^)' 


^/i  /'^+(d^H^V 


(2.46) 


it  is  evident  that 


n 


^^-^(aV+bs)  .  j  ^  j__(c,dJ,H^,d„,r) 


=  N  E  (^)”h  (— ^) 

/2  "  2/2a  (2.47) 


and  therefore  each  of  the  functions  can  be  found  where 


J„(o,dJ,H  d D  =  ^(^)“h  (-^) 
“  -LOO  /HT  /2  "  2/2a 


(2.48) 


And  finally,  these  functions  are  related  to  the  original 
integral  as 

/dH(Jjj(H)L(c,H)  =  S  f j  Jjj(c,d^,H^,d^,r) 

J 

In  order  to  find  the  functions  without  evaluating 
each  of  the  Hermites  directly,  let  hjjCX)  =  a**  H^(X)  so 
that  Eq  (2.42)  becomes 


2Xh^(X) 


.  2n  h^.iCX) 

a“^l  a“-l 


(2.49a) 
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or 

2(Xa)h^(X)  =  h^^3_(X)  +  2(na2)h^_^j^(X)  (2.49b) 

Now  if  ^  ~  (2.49b) 

is 

2(Xa)j^/HT  =  j^^^/(n+l)!  +  2(na2) j^_^,/(n-l)  !  (2.50) 

which  reduces  to 


j 


n+1 


avii  j 


n-1 


1 


(2.51) 


Since  the  functions  J  are  related  to  j  as 

n  *^n 


J^(c,dJ,H^,d^,r)  = 


N 


(/2) 


n  "n 


(2.52) 


can  be  redefined  as  J  '  =  Nj  and  with  H  =  1  and 
n  "n  n  o 

=  2X  ,  then  Jq'  =  N  and  =  2aXN  ,  so  that  the 

full  function  can  be  generated  from  Eq  (2.51).  Returning 
to  the  integral,  it  is  evident  that 


/dHO  (H)L(C,H)  =  I  f  .  — ii— 

j  ^  (/2)° 


(2.53) 


which  is  defined  as  r„(c). 

n 


^  Jn'  =  /dHffl  (H)L(C.H)  (2.54) 

°  (/2)°  j  J  “  “ 


Replacing  the  integrals  in  the  matrix  coefficients  of 
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Eq  (2.20)  gives 


B. 


nn 


P 


n 


q 


y 


(2.55a) 

(2.55b) 

(2.55c) 

(2.55d) 

Z 

(2.55e) 

Z 

(2.55f) 

It  is  evident  that  with  r„(c)  known,  all  of  the 

n 

matrix  coefficients  can  be  determined  with  the  exception  of 

and  C„„ .  The  smoothing  function 
nn  nq  ° 

/dH9)j^(H)Ag(-A2)9)^.(H')  (2.56) 


must  be  computed.  When  the  functions  are  assumed  to  be 

2 

Hermite  polynomials,  the  operator  g(-A  )  must  be  chosen  to 

provide  a  measure  of  roughness  that  relates  to  the  form  of 

these  functions.  Since  the  Hermite  polynomials  are  Harmonic 

oscillator  functions,  the  measure  of  roughness  can  be 

1  2  d^ 

appropriately  defined  by  the  operator  g(i(Z  -  — s))  as 

dZ'^ 
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(2.57) 


2 

dZ 

The  use  of  this  factor  in  the  solution  of  the  Harmonic 
oscillator  is  well  known,  where 

2 

|(Z^-  (2.58) 

Therefore,  the  roughness  term  is 

2 

=  A/dHdH-9>n(H)g(|(z2-  -^))0n^(H')  (2.59) 

dZ 

which,  upon  integration,  is 

"o  ’  “o«nn'8<”4)  -  “o«„n'  ”  <2-60) 

r*! 

By  defining  the  operator  as  shown,  the  solution  of  the 
smoothing  function  is  reduced  to  an  arbitrary  series  in 
(n+^)  consisting  of  R  terms.  The  svun  of  this  series  is 
added  to  the  diagonal  coefficients  of  The  quantities 

are  found  in  a  more  straightforward  manner  since  these 
coefficients  are  simply  the  Hermite  functions  evaluated  for 
the  arguments  (really  Hq/H^). 

Once  all  of  the  matrix  coefficients  are  known,  the 
matrix  problem  can  be  solved  for  the  "optimiim"  values  of 
Sq,  b^,  and  X^.  Using  these  values,  the  fit  of  the  magnetic 
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III.  Computer  Program  Development 

Two  separate  programs  were  written  according  to  the 
theory  discussed  in  Chapter  II.  These  programs  have  been 
designated  RPARAM  and  MOSFIT.  The  sole  purpose  of  RPARAM 
is  to  calculate  and  store  the  data  independent  integral 
values  called  in  Chapter  II.  Program  MOSFIT  solves 

the  minimization  problem,  then  calculates  and  plots  the 
magnetic  field  density  p(H).  Once  the  p(H)  is  appropriately 
modelled,  the  fit  spectrum  is  computed  and  plotted.  It  is 
evident  that  this  programming  can  be  divided  into  three 
steps : 

1.  Computation  of  1*^(0) 

2.  Determination  of  p(H) 

3.  Computation  of  Fit  Spectrum 

This  chapter  will  discuss  the  computer  program  development 
as  outlined  by  the  three  categories  above.  A  full  listing 
of  the  programming  may  be  found  in  Appendix  B. 

Computation  of  r^(c) 

Since  the  computation  of  the  parameters  is 

independent  of  the  actual  counting  data,  a  separate  program 
was  devised  to  calculate  and  store  the  values  of  r^^Cc)  for 
a  given  Mossbauer  nucleus.  The  flowchart  for  this  program 
is  shown  in  Figure  3,  Program  RPARAM  is  clearly  divided 
into  four  subroutines,  FIX,  CONS,  INTEG,  and  RITEO. 
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Subroutine  FIX.  The  program  RPARAM  requires  the 

input  of  the  values  ,  T,  f j ,  H^,  for  the  Mossbauer 

nucleus  studied.  These  values  are  material  dependent  and 

can  be  determined  by  examination  of  the  spectroscopic 

characteristics  of  the  nucleus  in  question.  The  values  for 
57 

Fe  have  been  computed  in  Appendix  A  and  are  shown  in 
Table  I. 


TABLE  I 
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Input  Parameters  for  Fe 


j 

d^(nm) 

f  r(mm/s) 

d^(mm/s) 

H^CkOe) 

1 

-.2437 

3  .23 

.258 

400 

2 

-.1413 

2  .23 

♦ 

* 

3 

-.0393 

1  .23 

★ 

4 

.0393 

1  .23 

* 

* 

5 

.1413 

2  .23 

* 

♦ 

6 

.2437 

3  .23 

♦ 

♦ 

♦  d 

It 

Q  and  are  independent  of  j 

is  clear  that  this  data  is  found 

in  a  variety 

of  units. 

The 

units  shown 

in  Table  I  are  those 

generally  found  in  the 

literature  for 

each  parameter.  For 

this  data  to 

be  useful 

in 

relation  to 

the  theory  developed 

in  Chapter  II 

,  all  of  the 

parameters  must  be  found  in  terms  of  channel  number.  The 
purpose  of  subroutine  FIX  is  to  convert  the  input  data  to 
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channel  data.  The  values  of  d^,  d^,  and  r  can  be  converted 
to  channel  data  with  two  conversion  factors.  Since  two  of 
the  parameters  are  already  in  units  of  millimeters/second, 
it  is  advantageous  to  first  convert  d^^  from  nuclear  magne¬ 
tons  to  mm/s.  This  can  be  accomplished  using  the  conver¬ 
sion  factor  shown 


CONV  =  3.1574417  E-12  KeV/nm-kOe  (3.1) 

in  combination  with  the  energy  of  the  nuclear  transition 
and  the  speed  of  light.  Therefore, 


d^(mm/s-kOe)  =  d^(nm) •CONV-c(mm/s)/E^(KeV)  (3.2) 


so  that  dj^H^  will  be  in  units  of  mm/s.  Finally,  each  of 
the  parameters  is  changed  to  channel  data  by  a  factor  that 
relates  velocity  to  channel  number.  This  factor  results 
from  a  least-squares  fit  of  channel  number  to  velocity 
developed  by  Captain  B.  E.  Pate  (Ref  13).  Captain  Pate 
found  a  given  velocity  to  be  related  to  channel  ntunber  as 
channel  =  (velocity  +  9.795990)/ .04819964 

where 

one  channel  =  .04819964  mm/s 
These  relations  are  incorporated  into  subroutine  FIX  to 
convert  from  velocity  to  channel  data. 

Subroutine  CONS.  A  determination  of  the  parameters 
r^(c)  first  requires  the  computation  of  several  terms  that 
are  used  in  this  calculation.  This  subroutine  solves  for 
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these  terms  as  outlined  by  the  flowchart  included  in 
Appendix  B. 

The  first  set  of  values  calculated  are  the  a  terms 
found  in  the  recursion  relation  (Eq  (2.51))  and  defined  in 
Eq  (2.36a).  Since  these  a  terms  vary  with  d^,  six  values 
are  computed  to  correspond  to  the  six  absorption  peaks  in 
the  spectrum.  The  evaluation  of  a  is  found  in  line  23  of 
RPARAM  as 

A(J)  =  SQRT(((D1(J)*H0)*#2 

-  GMMA(J)**2)/((D1(J)*H0)**2 

+  GMMA(J)**2))  (3.3) 

The  other  terms  to  be  evaluated  are  the  argument,  X,  of  the 

function  in  the  recursion  relation  as  well  as  the  con- 
“n 

stant  N.  The  argument  X  is  defined  in  Eq  (2.44) 

^  ^  diH„(d„-c) 


and  the  value  of  N  can  be  found  by  Eq  (2.46),  where 


N 


H^e  ^  r^+CdjH^)^ 

Vif  /’^+(d^H^)^ 


It  is  evident  that  each  of  these  terms  vary  with  both  the 
channel  number  and  the  absorption  peak  in  question.  The  X 
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is  found  in  line  27  of  RPARAM  by 

Z(J,CI)  =  D1(J)*H0*(D0-CI)/SQRT((D1(J)*H0)**4 

-GMMA(J)**4)  (3.4) 

The  constant  N  is  computed  in  two  parts.  NX  represents  the 
exponential  part  of  N  while  NC  combines  the  terms  not  asso¬ 
ciated  with  the  exponent.  As  a  result,  N  is  found  in  three 
steps: 


NX  =  EXP(-(CI-D0)**2/(2.*(GMMA(J)»*2 

+  (D1(J)*H0)**2))  (3.5a) 

NC  =  H0/(PI**(l./4. )*SQRT(GMMA(J)**2 

+  Dl(J)**2*HO**2))  (3.5b) 

N(J,CI)  =  NC*NX  (3.5c) 

Each  value  of  A,  Z,  and  N  are  used  in  turn  to  compute  the 
functions  and  the  r^(c)  parameters  in  subroutine  INTEG. 

Subroutine  INTEG.  This  subroutine  finds  the  func¬ 
tions  by  means  of  the  recursion  relation  presented  in 
Eq  (2.51),  where  Jq  and  are  known.  Each  of  these  func¬ 
tions  will  vary  with  channel  number  and  absorption  peak. 

The  recursion  relation  is  found  in  line  47  of  RPARAM  and  is 
coded  as 

FNJ(J,S+1,CI)  =  (2.*A(J)/SQRT(S+1))*(J,CI)*FNJ(J,S,CI) 
-  A(J)*FNJ(J,S-1,CI)*SQRT(S1))  (3.6) 
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The  values  of  r^(c)  are  subsequently  determined  in  line  50 
where 

R(S,CI)  =  F(J)*(1./((SQRT(2. ))**S)*FNJ(J,S,CI) 

+  R(S,CI)  (3.7) 

which  corresponds  to  Eq  (2.54).  A  flowchart  of  this  routine 
can  be  found  in  Appendix  B. 

Subroutine  RITEO.  The  function  of  RITEO  is  simply  to 
write  the  values  of  r^(c)  for  permanent  files  storage.  It 
is  designed  to  write  a  maximum  of  five  values  per  line  in 
an  E14 . 7  format . 

Determination  of  p(H) 

The  magentic  field  density  is  determined  in  MOSFIT 
using  values  found  in  the  least-squares  minimization.  The 
subroutines  of  interest  in  this  case  are  READN,  COEFF, 

MATRl,  HERMIT,  MESCHP,  MAGFD  and  PLTMAG.  A  generalized 
flowchart  of  the  entire  MOSFIT  program  is  found  in  Figures 
4  and  5. 

Subroutine  READN .  Subroutine  READN  is  designed  to 

read  the  experimental  data  and  the  previously  computed 

values  of  r^(c)  from  a  data  tape.  The  data  is  read  from 

TAPED  in  an  I3,10F7,0  format.  This  format  corresponds  to 

the  data  as  recorded  by  the  multichannel  analyzer .  The 

r„(c)  values  are  read  from  TAPES  in  format  chosen  to  corre- 
n 

spond  with  RITEO. 
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Subroutine  MATRl .  The  utility  of  this  routine  is  to 
assign  the  previously  computed  coefficients  to  a  matrix  A 
so  the  matrix  problem  AX  =  Y  can  be  solved.  This  routine 
could  be  incorporated  into  the  subroutine  COEFF  if  desired. 

Subroutine  HERMIT.  Subroutine  HERMIT  calculates 
and  stores  the  n+1  Hermite  polynomials  used  in  the  compu¬ 
tation  of  p(H).  The  range  of  magnetic  field  strength  is 
divided  into  100  sub-divisions  and  the  set  of  polynomials  is 
evaluated  for  each  increment  of  H/H^.  The  polynomials  are 
generated  using  the  known  values  of  H^  and  H^^  and  the 
Hermite  recursion  relation  in  Eq  C2.42) 

2XH^CX)  =  +  2nH^_3^(X) 

This  relation  is  coded  in  line  174  of  MOSFIT  as 

HPCI,S+1)  =  2.*CHCI)/H0)*HPCI,S) 

-  2.*S*HPCI,S-1)  (3.10) 

where  the  index  I  represents  the  increment  of  the  field 
strength  and  S  represents  the  polynomial  number  from  0  to  n. 

Subroutine  MESCHP .  Subroutine  MESCHP  solves  the 
matrix  problem  by  performing  a  Choleski  decomposition  on 
matrix  A.  (Ref  14).  The  program  was  written  by 
Dr.  Shankland  and  is  found  in  a  library  cataloged  as  BINDEK, 
CY=1, ID=T760142  in  the  AFIT  system.  The  arguments  for  the 
subroutine  are: 
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A  -  the  coefficient  matrix  of  the  problem 
NA  -  the  dimensioned  size  of  matrix  A 
K  -  the  actual  size  of  A  (KiNA) 
y  -  the  solution  of  matrix  Y 

1  -  indicates  Y  is  a  vector  dimensioned  CNAxl) 

IE  -  error  indicator  (no  input  needed) 

D  -  scratch  vector  must  be  dimensioned  to  at  least 
length  K 

The  answers  to  the  matrix  problem  are  found  in  Y  at  the 
completion  of  the  computation. 

Subroutine  MAGFD .  This  routine  uses  the  solution 
vector  and  the  Hermite  polynomials  HP{I,S)  to  model  the 
magnetic  field  density.  Since 

pCZ)  =  S  H^(Z) 
n 

and  the  expansion  coefficients  a^  are  found  in  the  solution 
vector  Y,  it  follows  that 

RHO(I)  =  YCS)»HPCI,S)  +  RHO(I)  (3.11) 

when  the  equation  is  summed  over  S,  Equation  (3.11)  is 
found  in  line  197  of  MOSFIT  with  the  addition  of  constants 
for  orthogonality  and  normality.  This  equation  is  used  to 
find  the  value  of  p(Z)  at  the  incremented  value  of  Z. 

Subroutine  PLTMAG .  Simply  put,  PLTMAG  uses  library 
software  (the  DISSPLA  package)  to  plot  the  values  of  RHOCI) 
versus  the  incremented  values  of  the  magnetic  field  strength 
HCI).  It  requires  DISSPLA , ID=LIBRARY, SN=ASD  be  attached  as 
a  library. 
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The  solution  of  the  magnetic  field  density  can  be 
more  clearly  understood  by  examining  the  flowcharts  for 
COEFF,  HERMIT,  and  MAGFD  that  are  contained  in  Appendix  B. 

Computation  of  Fit  Spectrum 

With  the  bulk  of  the  calculations  done,  the  final 
step  is  to  compute  the  fit  spectrum.  The  fit  spectrum  is 
found  using  the  RHO(I)  previously  computed  and  the  single 
subroutine  SPCTRM. 

Subroutine  SPCTRM .  The  fit  spectrum  can  be  found  in 
two  ways.  Once  the  minimized  parameters  a^^,  b^  and  have 
been  found  in  the  solution  vector  Y,  these  values  can  be 
applied  in  Eq  (2.2)  to  find  the  fit  spectrum.  The  integral 
is  approximated  by  the  rectangular  method: 

AH  Z  p(H«)L(c,H«) 

OC 

using  the  p(H)  calculated  in  MAGFD  and  the  lineshape  formula 
found  in  Eq  (2.24).  The  results  is  coded  in  MOSFIT  as 

RL(I)  =  F(K)*EXP(EX)/(GMMA(J)*SQRT(2.*PI))*DELH*RH0CJ) 

+  RL(I)  (3.12) 

where  each  of  the  parameters  r,  d^^,  and  d^  are  corrected 
by  subroutine  FIX  as  in  RPARAM.  An  alternate  method  for 
determining  the  spectrum  has  also  been  programmed  where  the 
integral  values  rjj(c)  and  the  computed  coefficients  a^  are 
used.  Since 
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/p^(H)L(c.H)dH  =  T^(c) 

the  integral  part  of  Eq  (2.2)  can  be  found  as 

^  Vn<‘=> 

n 

so  that  the  parameters  T,  d^,  d^,  and  subroutine  FIX  are 
not  necessary.  This  method  has  two  advantages: 

1.  Computation  time  is  greatly  reduced. 

2.  The  analytical  solution  in  r^(c)  is  expected  to 
be  more  accurate  than  the  corresponding  numerical  integra¬ 
tion. 

An  outline  of  the  fitting  method  using  the  r^(c)  values  is 
flowcharted  in  Appendix  B. 

Subroutine  PLTSPT.  This  routine  also  uses  the 
DISSPLA  package  to  plot  the  fit  spectrum.  The  fit  points 
RL(I)  are  plotted  versus  the  channel  number  C.  For  compara¬ 
tive  purposes,  the  actual  data  N(C)  is  printed  on  the  same 
plot . 

Other  Routines 

In  addition  to  the  programming  previously  discussed, 
there  are  two  essential  subroutines  that  have  not  been 
written  but  are  necessary  for  the  solution  of  this  problem. 

Subroutine  CONST.  This  subroutine  is  necessary  for 
the  application  of  constraints  in  the  problem.  The  principle 
purpose  of  CONST  would  be  to  calcdlate  the  matrix  coefficients 
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C  found  in  Eq  (2.55).  The  computed  points  RHO(I)  are 
analyzed  for  negative  values  and  the  index  I  of  the  points 
to  be  constrained  are  noted.  Using  the  Hermite  values 
already  computed,  the  coefficients  C  are  simply  the  n 
Hermite  values  of  each  negative  point 

Subroutine  MATRN .  MATRN  would  be  analogous  to  sub¬ 
routine  MATRl;  however,  it  would  only  assign  the  coeffi¬ 
cients  C  to  the  matrix  A.  Matrix  is  dimensioned  to  allow 
the  addition  of  constraints  as  necessary. 

Besides  subroutines  CONST  and  MATRN,  an  additional 
routine,  SMOOTH,  has  been  written.  SMOOTH  functions  as  a 
test  of  accuracy  in  the  fitting  procedure.  This  routine 
calculates  the  value  of  the  system  S  (Eq  (2.12))  after  the 
least-squares  minimization.  This  value  is  used  to  determine 
when  a  sufficient  number  of  Hermit  polynomials  has  been  used 
to  describe  p(H).  This  procedure  is  detailed  in  Chapter  V. 
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IV.  Program  Validation 

In  order  for  a  computer  program  to  effectively  solve 
a  desired  problem,  the  program  must  be  examined  to  determine 
whether  it  is  functionally  and  computationally  accurate.  An 
excellent  discussion  on  the  framework  and  procedures  of  pro¬ 
gram  testing  can  be  found  in  the  literature  (Ref  15). 

The  individual  subroutines  of  both  programs  were 
examined  to  determine  the  correctness  and  accuracy  of  the 
numerical  calculations.  In  addition,  the  function  of  non- 
computational  routines  was  investigated  to  assure  desired 
performance.  The  current  programming  is  a  result  of  correc¬ 
tions  made  through  this  somewhat  iterative  testing  procedure. 
Therefore,  corrected  errors  will  not  be  discussed.  However, 
validation  procedures  and  numerical  testing  of  these  pro¬ 
grams  are  included. 

Program  RPARAM  was  initially  examined  with  regard  to 

the  input  data  (as  described  in  Appendix  A)  and  the  line- 

shape  function.  It  was  important  for  the  data  and  the 

lineshape  function  (Eq  (2.21))  to  produce  a  curve  comparable 
57 

to  an  Fe  spectrum.  This  test  also  helps  support  the 
correction  factors  used  in  subroutine  FIX.  The  result  of 
this  examination  of  data  and  lineshape  can  be  seen  in 
Figure  6.  This  plot  was  generated  using  the  data  in  Table  I 
(as  corrected  by  subroutine  FIX).  The  general  lineshape 
correlates  with  the  expected  results. 
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with  Given  Data 


Next,  the  computational  accuracy  of  RPARAM  was 
examined.  It  was  arbitrarily  decided  to  examine  the  values 
of  all  parameters  evaluated  at  J=l,  CI=200  and  S=1  to  4. 

The  relations  FNJ  (corresponding  to  were  hand  calculated 
by  the  following  two  equations; 

A(J)®H  (Z) 

fnj  = - ^ - 

s 

(a) 


FNJ 


s+1 


=  2A(J) 
/s+1 


[Z(FNJg)  -  A/E  (FNJg_^)] 


(b) 


and  computationally  determined  by  the  programmed  recursion 
method  as  well  as  from  the  actual  series  (analytic)  solution 
(Eq  (2.40)).  The  results  were  vaforable  as  shown  in 
Table  II. 


TABLE  II 


Evaluation  of  FNJ 


s 

FNJ  (a) 

FNJ  (b) 

Computer  Recursion 

1 

-.1298061 

-.1298061 

-0.1298062505182 

2 

-1.3986532 

-1.3986533 

-1.39865339369 

3 

.3162459 

.3162459 

0.3162465049902 

4 

2.395768 

2.3957681 

2.3957683649 

In  addition,  A(l)  was  found  to  correspond  to  the  calculated 
value  of  .9987102  while  Z  varied  slightly  in  the  third 
decimal  place. 
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The  final  subroutine  RITEO  was  assumed  to  be  func¬ 
tionally  correct  based  on  examination  of  the  data  tape 
generated. 

Program  MOSFIT  was  also  examined  for  proper  perfor¬ 
mance.  Subroutine  READN  was  tested  simply  by  writing  out 
the  data  read  into  the  program  and  comparing  the  printed 
results  with  the  actual  data  files.  Subroutine  COEFF  was 
tested  with  simple  integer  values  for  r^^Cc)  and  10  channels. 
Each  matrix  coefficient  was  calculated  by  hand  and  compared 
to  the  results  from  COEFF  using  the  same  data.  The  computed 
values  were  found  to  be  in  agreement  with  those  that  were 
hand  calculated.  The  routine  MATRl  was  tested  by  printing 
out  the  values  of  A  and  Y  in  the  matrix  problem  AX  =  Y 
The  matrix  A  was  symmetric  as  expected  and  the  coefficients 
corresponded  correctly  with  those  computed  in  COEFF. 

HERMIT  was  found  to  correctly  compute  the  Hermite  polynomials 
H^(H/H^)  based  on  a  comparison  between  the  computed  Hermites 
and  those  hand  calculated  from  the  known  polynomial  forms, 
i  .e. , 

H^CZ)  =  1 

H^(Z)  =  2Z 

HgCZ)  =  4Z^  -  2 

Hg(Z)  =  8Z^  -  12Z 

etc. 

One  of  the  fundamental  subroutines  is  MESCHP  which  solves 
the  matrix  problem  using  a  Choleski  decomposition  of 
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matrix  A.  This  routine  was  tested  on  a  simple  2x2  matrix 
problem  to  ascertain  its  computational  accuracy.  Using  the 
problem 


2 

1 


the  routine  solved  x  as  1  and  y  as  2,  which  is  correct. 

The  subroutine  MAGFD  was  analyzed  by  using  a  single 

Hermite  to  model  the  magnetic  field  density  so  that  the 

plotted  p(H)  should  correspond  to  the  known  shape  of  the 

-z^  /2 

Hermite  used.  Using  and  the  orthogonality  factor  e  '  , 
the  subroutine  generated  a  Gaussian  centered  at  0,  corre¬ 
sponding  to  a  plot  of  the  Hermite.  Neither  plot  routine  was 
tested  independently,  but  in  the  course  of  other  test  pro¬ 
cedures,  the  plotting  functioned  successfully. 

Finally,  the  functional  accuracy  of  the  entire  pro¬ 
gram  was  tested  using  data  with  a  known  p(H).  This  data 
was  generated  using  a  magnetic  field  density  of 


p(Z)  =  r  e~^^/^(H^(Z)  +  I  HgCZ)) 


and  a  linear  background  of  100  counts  per  channel.  The 
resulting  spectrum  and  magnetic  field  density  can  be  seen  in 
Figures  7  and  8,  respectively.  The  test  procedure  used  this 
generated  counting  data  and  the  previously  computed  r^(c) 
values  in  program  MOSFIT.  It  is  evident  that  if  the  pro¬ 
gram  is  operating  properly,  the  plots  resulting  from  the 
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least-squares  minimization  in  MOSFIT  should  be  virtually 
identical  to  those  in  Figures  7  and  8.  However,  as  evidenced 
by  Figures  9  and  10,  the  plots  do  not  correspond.  The 
expansion  coefficients  for  the  magnetic  field  in  Figure  9 
were  found  to  be 

aQ  =  0.4842 

ag  =  0.0864 

which  are  roughly  one-half  the  value  of  the  input  coeffi¬ 
cients  aQ  =  1  and  ag  =  0.1667  .  Although  the  background 

count  of  ~102  counts  does  roughly  correspond  to  the  100 
counts  of  the  input,  it  appears  that  a  factor  of  two  is 
missing  somewhere  in  the  program.  Therefore,  the  expansion 
coefficients  were  multiplied  by  two  and  the  program  was  run 
again.  The  results  are  shown  in  Figures  11  and  12.  These 
results  reproduce  the  general  shape  of  p(H)  and  the  fit 
spectrum  seems  to  match  the  test  input  rather  well.  In 
addition,  the  fit  spectrum  was  determined  using  the  r^(c) 
values  rather  than  a  numerical  integration.  An  apparent 
discrepancy  occurs  in  these  runs.  The  fit  spectrum  found 
without  multiplying  the  expansion  coefficients  by  two 
follows  the  actual  data  more  closely  than  the  "corrected" 
spectrvim.  This  can  be  seen  in  Figures  13  and  14.  The  mag¬ 
netic  field  densities  corresponding  to  these  results  are 
identical  to  those  in  Figures  9  and  11.  Therefore,  it  is 
obvious  that  although  the  spectrum  in  Figure  13  appears  to 
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Figure  9.  Magnetic  Field  Density  from  Initial 
Test  Run 
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in  Initial  Run 


RHO(H)  FOUND  FROM  GENERRTED  DflTfl 
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Figure  11.  Magnetic  Field  Density  With 


FIT  SPECTRUM  VS.  RCTUflL  DRTfl 


and  2x  Factor 


SPECTRUM  USING  R  VRLUES-UNCORRECTED 


0'S9 

SINHOO 


o-sv  o'se 


Figure  13.  Fit  Spectrum  Using  r  (c)  Values 
Without  2x  Factor 


riT  SPECTRUM  USING  R  VRLUES-CORRECTED 
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Figure  14.  Fit  Spectrum  Using  Values 

and  2x  Factor 


be  a  "good"  fit,  the  corresponding  magnetic  field  density 
does  not  match  the  input.  These  problems  were  not  resolved. 
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V.  Results ,  Recommendations  and 
Conclusions 

The  objectives  of  this  study  have  not  been  entirely 
met  because  of  an  unknown  error(s)  in  the  computer  imple¬ 
mentation  of  this  theory.  This  chapter  will  examine  the 
proposed  optimization  of  the  results  in  a  fully  operational 
program,  conclusions  drawn  from  the  work  done  to  this  point, 
and  recommendations  for  further  study  and  modification. 

Optimizing  Results 

It  is  obvious  from  Chapter  IV  that  no  reasonable 
results  have  been  produced  since  no  actual  data  has  been 
processed.  However,  in  the  fitting  of  p(H),  the  question 
arises  as  to  the  niunber  of  orthonormal  functions  required 
to  "best-fit"  the  system.  This  problem  is  solved  if  the 
system  S  (Eq  (2.12))  is  evaluated  after  the  minimization. 

That  is,  using  the  minimized  parameters,  a^^,  b^ ,  X^,  find 
a  value  of  S  for  the  given  number  of  polynomials  n.  The 
expected  result  of  a  plot  of  InCs)  versus  n  is  shown  in 
Figure  15,  This  figure  shows  that  the  value  of  S  is  expected 
to  decrease  with  increasing  n,  then  level  off  in  a  plateau 
region,  followed  by  a  sharp  drop  off  at  large  n.  The  best 
fit  is  found  with  an  n  in  the  plateau  region  where  S  is  at 
its  minimum  value.  The  region  of  sharp  decrease  in  S  after 
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the  plateau  is  not  useful  because  the  fit  begins  to  follow 
the  actual  data  too  closely  beyond  this  point. 


ln(S) 


Figure  15.  Result  Optimization 


Conclusions 

The  problem  with  the  programming  is  not  known; 

however,  an  error  of  a  factor  of  two  should  be  able  to  be 

found  in  a  straightforward  manner.  First,  in  considering 

the  actual  input  data  (Figures  7  and  8),  it  was  noted  that 

the  magnetic  field  density  was  calculated  without  the 

1  2 

normalization  constant  ( - )  for  H  .  This  constant  has 

a  value  of  approximately  0.75.  Since  the  program  MOSFIT 
uses  the  orthonormal  functions  55^  in  the  expansion  of  p(H), 
this  is  clearly  one  source  of  disagreement  between  input 
and  results.  The  discrepancy  that  arises  when  the  r^(c) 
values  are  used  is  more  difficult  to  explain.  These  r 
values  are  computed  with  a  factor  of  (1 .  / .  This  factor 
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reduces  the  magnitude  of  the  second  order  r  terms  by  one- 
half.  Obviously,  if  this  factor  is  incorrect,  another 
source  of  disagreement  is  found. 

An  interesting  observation  that  has  been  made  is  that 
all  odd  r  terras  are  zero  and  can  be  ignored  in  these  calcu¬ 
lations.  This  fact  is  intuitively  obvious,  based  on  the 
early  assumption  that  the  magnetic  field  density  is  symmet¬ 
ric  about  zero.  In  addition,  with  careful  examination  of 
Eq  (2.35b)  and  the  definition  of  the  Hermit  argument  in 
Eq  (2.43),  it  is  easily  shown  that  all  ^^^(c)  values  must 
equal  zero  for  odd  n. 

It  is  also  significant  to  note  that  although  Pate's 
velocity  fit  is  dependent  on  specific  data,  it  does  not 
necessarily  eliminate  the  data  independent  character  of  the 
rn(c)  values.  The  velocity  fit  corresponds  to  an  NFE 
spectrum,  which  is  used  as  a  "base  spectrum"  for  the  deter¬ 
mination  of  the  magnetic  field  density.  Since  this  same 
base  spectrum  will  be  used  each  time,  the  velocity  param¬ 
eters  do  in  fact  become  data  independent.  However,  if  a 
specific  velocity/channel  fit  is  desired  for  each  set  of 
data,  this  can  be  readily  incorporated  into  the  program. 

The  multichannel  analyzer  that  stores  the  counting  data 
includes  velocity  data  as  well. 
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Finally,  these  results,  although  far  from  conclusive, 
appear  to  show  that  with  minor  corrections  this  program  may 
provide  a  viable  method  to  fit  the  magnetic  field  distri¬ 
bution  from  amorphous  metal  alloys. 

Recommendations 

The  principal  recommendation  must  be  that  this  study 
be  continued  to  conclusion  with  other  parties  investigating 
the  theory  and  programming  in  order  to  find  the  minor  error. 
Simple  manipulation  of  the  existing  programs  should  ascertain 
if  the  two  factors  discussed  earlier  are  contributors  to  the 
error.  Obviously,  the  current  programming  should  be  modi¬ 
fied  to  eliminate  the  computation  and  use  of  all  odd  r^(c) 
values,  since  all  of  these  values  are  zero.  In  addition, 
an  important  change  would  be  to  find  a  better  way  to  do 
channel  to  velocity  conversion.  If  the  entire  programming 
could  be  redeveloped  in  relation  to  velocity  rather  than 
channel  number,  conversion  routines  would  not  be  necessary 
and  computation  time  would  be  enhanced.  Also,  in  determining 
the  fit  spectrum,  it  is  recommended  that  the  r^(c)  parameters 
be  used  instead  of  numerical  integration  routines.  This  will 
increase  the  speed  and  should  improve  the  accuracy  of  the 
results. 

Finally,  any  computational  problems  in  the  solution 
of  the  matrix  problem  may  be  overcome  by  using  double  pre¬ 
cision  numbers.  However,  with  the  large  storage  space 


already  in  use,  the  longer  turn-around  time  may  negate  any 
benefit  gained.  Another  interesting  facet  would  be  to  add 
a  statistical  Goodness  of  Fit  test  to  the  programming  in 
order  to  test  the  final  spectrum.  Finally,  once  this  pro¬ 
gram  is  fully  operational,  it  would  be  interesting  to  examine 
its  effectiveness  on  spectra  of  amorphous  metal  alloys 
since  the  programming  was  developed  for  this  purpose  and 
these  materials  have  been  of  recent  interest  to  the  Air 


Force. 
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Appendix  A 

Input  Parameters  for  Program  RPARAM 


Program  RPARAM  requires  the  input  of  certain  param¬ 
eters  that  describe  the  expected  spectral  shape.  For 
example,  note  the  relationship  between  these  parameters  and 
the  spectrum  in  Figure  A-1. 


Figure  A-1.  Relation  of  Input  Parameters 
to  Spectrum 


The  purpose  of  this  discussion  is  to  direct  the  reader  to 

the  correct  quantities  in  the  literature  for  input  into 

this  program.  Since  the  absorber  in  amorphous  metal  alloys 
57 

is  Fe  ,  the  quantities  related  to  this  isotope  will  be 
discussed  here. 

Two  of  the  quantities,  T  and  are  of  a  rather 
arbitrary  nature.  Gamma  (T),  which  represents  the  spectral 
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linewidth  can  be  taken  as  a  quantity  somewhat  larger  than 
the  natural  linewidth  of  the  nucleus  under  investigation. 
This  follows  from  the  expected  broadening  of  the  peaks  with 
respect  to  the  natural  linewidth.  The  linewidth  was  taken 
as  .23  mm/sec  based  on  previous  work.  The  natural  line- 
width  is  found  as  .194  mm/sec  in  the  literature  (Ref  2). 

is  an  arbitrary  scaling  factor  for  the  dimensionless 
variable  Z.  is  chosen  to  give  an  appropriate  distribu¬ 

tion  p(Z).  If  is  chosen  too  large,  the  distribution  is 
spread  over  a  wide  range  of  field  strength,  while  a  smaller 
will  narrow  the  distribution.  The  in  this  program 
was  varied,  but  the  best  results  seem  to  be  found  with 
=  400  kOe. 

The  f  factor  is  not  the  recoil-free  fraction  found 
in  the  literature.  Instead,  f  represents  the  relative 
intensity  of  each  ..p'^ctral  peak.  The  calculated  intensity 
of  the  individual  components  of  a  spectrum  involves  rather 
complex  calculations  and  can  be  found  in  the  literature 

57 

(Ref  1).  The  relative  intensities  of  the  lines  in  the  Fe 
spectrum  integrated  over  all  directions  is  found  to  be 
3;2:1:1;2:3  (Ref  2). 

Finally,  the  spreading  quantities  are  best  understood 

57 

by  examining  the  splitting  of  the  Fe  nucleus  (Figure  A-2). 
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Figure  A-2.  Splitting  of  the  Fe  Nucleus 
The  energy  eigenvalue  of  a  given  level  m  is  found  as 

where  gj^  is  the  nuclear  Lande  splitting  factor,  6j^  is  the 
nuclear  Bohr  magneton  and  mj.  is  the  magnetic  quantum  number. 
However,  the  nuclear  magnetic  moment  is 

V  =  CA-2) 

and  the  energy  eigenvalues  are 

-ym^ 

=  C-^)H  CA~3) 

which  are  evaluated  in  Figure  A-2.  The  AE„  for  a  transition 

m 

is  found  from  the  difference  in  E„  for  the  allowed  transi- 

m 

tions.  Therefore,  the  energy  of  a  gamma-ray  is  found  as 
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(A-4) 


The  energy  is  perturbed  by  the  isomer  shift  (I.S.)  so 
that 

E^  =  +  I.S.  +  (A-5) 

Y  o  m 

From  Eq  (A-5),  it  follows  that  and  d^^  are  found  as 

d„  =  E  +  I.S. 
o  o 

di  =  AE„/H 
1  m' 

Since  E^  relates  to  a  zero  velocity,  d^  =  I.S.  ,  which 

gives  the  input  to  RPARAM,  when  the  values  of  and  are 
57 

known.  For  Fe  >  . 153nm  and  =  +.0903nm 

It  is  evident  that  there  is  flexibility  in  the  defi¬ 
nition  of  the  d  parameters  discussed  above.  This  study  has 
only  considered  the  magnetic  field  interactions  in  the 
Mossbauer  analysis.  If  electric  field  or  quadripole  inter¬ 
actions  are  of  importance,  these  effects  could  be  included 
in  the  d  terms. 
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Appendix  B 

Computer  Program  Listings 

The  computer  programming  for  this  problem  was  written 
in  FORTRAN  V  for  a  CDC  Cyber  mainframe.  Programs  RPARAM 
and  MOSFIT  are  discussed  in  detail  in  Chapter  III  and  will 
not  be  examined  here.  However,  it  should  be  noted  that 
MOSFIT  does  require  the  addition  of  one  card  for  input  data. 
The  program  determines  constraints  based  on  user  input. 

All  other  data  is  entered  in  DATA  statements.  Also, 
plotting  requires  the  DISSPLA  package  and  online  plotting 
capability.  Program  DATA  is  a  modification  of  subroutine 
SPCTRM  in  MOSFIT  and  was  used  to  generate  data  for  program 
testing . 

Flowcharts  of  selected  subroutines  are  also  included 
in  this  Appendix.  Subroutines  CONS,  INTEG,  COEFF,  HERMIT, 
MAGFD  and  SPCTRM  have  been  examined.  The  flowcharts  for 
these  routines  wi^l  be  found  immediately  following  their 
individual  listings. 

The  program  listings  that  follow  are  numbered  for 
ease  in  finding  the  subroutines.  The  index  is  as  follows: 

RPARAM 


Subroutine 

Page 

Line  Number 

CONS 

73 

17 

INTEG 

75 

33 

RITEO 

77 

56 

FIX 

78 

72 

70 


MOSFIT 


Subroutine 


READN 

COEFF 

MATRl 

HERMIT 

MAGFD 

SMOOTH 

PLTMAG 

SPCTRM 

PLTSPT 


Page 

Line  Number 

80 

39 

81 

61 

85 

121 

86 

154 

88 

180 

90 

203 

91 

239 

92 

260 

96 

304 
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381  = 

RETURN 

Z82= 

ENC 
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Figure  B-6.  Flowchart  tor  Subroutine  MAGFD 
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li:-- 

£.£55'J7!V: 

'ai-. 

IN*£C.5  3i3I 

'!iC. 

Lit-.  - 

CG''’<;'./lNL7/K.L.tl 

Z'2o- 

CC'”'.:\/’N?T/5.g:3.488).N(;88) 

Z07-- 

DIMSl-lIOfii  rtf.4).£f£t£) 

LV^J“ 

DSTS  A/1./ 

£89=: 

SU«2=8. 

218= 

DO  18  CI=1>K 

«  . 

L*  .  “ 

SU!'1=8. 

21Z= 

DO  £8  S--8.L 

£’.3= 

s;jy,i=TisnR(S.ci»+S!j>ii 

£14=28 

CONTlfiOE 

e. 

3iJ.r£=S!i::)i(X!CI)-T(l+2!-r(L+3)*CI 

£li= 

!  «£+S'0«2 

£17=18 

CCNTINOE 

£15= 

CONS’sg, 

c:  :/  virito 

£28= 

ftl=L+4 

£21= 

£0  40  5=8fL 

H.L^ 

00X67=1  IKL+fiS)*!  ;S>4C(S.NQ)+C0N3T 

££?=48 

CO^tTIKUE 

£;i=38 

C’JNT  livjE 

tiC- 

ROOCH=8. 

££i= 

EG  58  6=8. L 

i.LI~ 

DO  60  66=8.1 

_ 

ti:'= 

IF  (S.NE.55)  GOTO  68 

^«o- 

tl.7- 

DO  78  1=1. IR 

4.01/* 

A4=EJP(-<!St.5!/(£.»It)n21 

'>•71  . 

RC=AA4(l+(S+.5)i«£H»I+R0 

£32=70 

COMIWJE 

RG.GH=«HT>'S>«2»R04RCyCH 

£34=i3 

CCNTIROE 

£35=58 

CCNTINliE 

236= 

FIT=,5*SUf!l-C0SST-rA*R0UCH 

237= 

RETURN 

£36= 

END 
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S:  PLT''iG!NS> 

THIS  SLSRC^TINE  UStS  THE  CIS5PLA  PACKAGE 
TO  P.OT  THE  mhVlt  rIE:.S  CESSITt  AS  A 
F'JNCTIOS  0*  K.  THE  CISSFLA  PKG  RUST  BE 
ATTACHED  AS  A  LIEEART  TO  USE  THIS  ROUTINE. 

CALL  CCRPRS 
CALL  SGNPLINSt 
CALL  GRACElg.l) 

CALL  TirttCSMCIH)  POUND  PRC.’'  CE‘«ERS'ED  CATAJ’,-‘e0i 

pie.:  H»Mgg.'Rr'C:H>»Mg2.7.5.4.5) 

CAL.  SRAPH<g..202..2...25) 

CALL  CL'VEIH.SHO.IM.g) 

CA..  :';i’L(Sfi) 

NKsNfi+1 

RETURN 

ENG 


f 
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S:tR:;T!'.i  iRC’Rs'tT.NA) 

CC^“C-c':SV£/F!6)  ,Ci  Ut  .SLi4?0l  .GKrA(4)  .D8 
EI!<EN?I3S  T;!iAJtEKGR«{400t 

2 -  a  C^?'.A.£2.Fi; .23,  .23, .23, .23, .23, .23, .24,3. Ii55?27/ 

SATA  F,Sl/3.,2.,l..l.,2.,3.,-.2437,-.ltl3,-.«393,.0393, 

:. 1413,. 2437/ 

DATA  '!\kth.t 

.y,.;  n-t'-'v:  %3  3"£CrRJ« 

“iv'.3':C  "rELO 

CtVilTf  UsI.‘.S  >Hi  2A■.'^C•^ .  '■’.-■-..ir. 


2S5-- 

284=27 


CAwL  r:)(IK'<,D0,Cl,C!'i'A) 

r:  18  i--i,« 

EXDRNS !  n '-r  tT  (143?*:-'' ;l+4h:m2 
F-;;i=2. 

DC  25  IKsl.KK 
13  27  3=1,108 

El={-(l-50-(Ol{!KUHyn)M2/(2.«C!<flAaK)«21) 

IF  (E7.LT.-250.)  THEN 
ftl(!)=RLlI)+0. 

ELSE 

RL!I)=F(IK>»E)(P(E)()/IC!!?!A(IK)4S3RT<2.i?:)i4DELKiW(3)4RLtI) 

EJiDIF 

C0STI«iL‘E 


2eS=  RLlDtEKGRNDdt-RLII) 

289=10  CONTINUE 

290=  RETURN 

291=  END 


// 
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Figure  B-7.  Flowchart  for  Subroutine  SPCTRM 
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:=:=  rDCK.'^si.cxnA) 

Z9:=  D:»ZNs:a-;  CKfo.G'-irAiKi 

2viz  li-'A  E?,::--SE:/:-.-llZ5.2.997925£n/ 

295=  DATA  CZ'.v.'S.lSTAi'JMZ/ 

29E=  CC  18  J--1.K 

297=  El  ( J>  =C0NV»D1 !  J)iC!(»;5'C/tS 

198=  Dl(j:=31!JI/.8i81996A 

299=  G>!''S(.Ji=C^^A!J)/.8;SI99t4 

188=18 

281=  E8=!E3+9.795998W.««199i4 

382=  RETURN 

303=  EKE 
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30M  SUBROUTINE  ?ITSPT(N?I) 

305*  REAL  S 

304*  CC»^CN/INPT/!!:0:3.A03!.S(M«) 

307*  COH?ICM/CRVE/r  14;  ,D1 14)  .RL(40'J5  .CMA(4)  .L3 


308= 

EIJ'ENSICSI  Cli2i) 

307= 

SO  20  K=1.40S 

Oi»- 

C1K)=K 

311=20 

CGriN'JE 

0  i  f 

CALL  EC'iPLCiN; 

313= 

CAL-  TI'LEi'R::  S^EC'RJ?!  VS 

314= 

■.:22.’CeUNT3IM30.7.5.t.5; 

p- 

V  4  V  “ 

CALL  CR177:P..55..i:..20.) 

314= 

CAv.  CL3VEi3.RL.«0.?; 

317= 

CALL  ff  ASKER  ill 

513= 

CA-L  •"LPIClJCi; 

319* 

I-  A 1.  i.  C  L  S  V  E I C  t  N 1 A  2  5 1  - 1 1 

320* 

CA-L  SESE'I’SCLPIC’l 

•J71- 

iLl- 

CALL  GRACE!. 01! 

^Li.“ 

CA.L  E'ACE 

CLw* 

CA-L 

vt4- 

CALi.  DCN£:l 

1.  wv* 

RETURN 

324= 

ENC 

AOTUSL  I;ATA»'»-1^0.’CHA>!-iEL 


96 


K6--  KK=0 
41=  DC-  40  1  =  1,40 

42=  lIRITEdS, 10001  KK,(RL(i:),n=K^+l.KK+l0) 

43=  KK=KK+14 

44=4f  CGN'INIE 

45=1000  FC-:»ATd:M3?7.0l 

H--lied  FCR'!4T(13,5'l4.fc> 

17=  ?TCP 

lf=  L>iB 
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